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Abstract 


We  have  analyzed  the  transferability  of  a  previously  proposed  Buckingham 
repulsion-dispersion  intermolecular  potential  for  the  explosive  hexahydro- 1 ,3,5-trinitro- 1 ,3,5-s- 
triazine  (RDX)  (Sorescu,  D.  C.,  B.  M.  Rice,  and  D.  L.  Thompson,  Journal  of  Physical 
Chemistry  B,  vol.  101,  p.  798, 1997)  to  predict  the  crystal  structures  (within  the  approximation 
of  rigid  molecules)  of  a  database  of  30  nitramines.  These  include  acyclic,  monocyclic,  and 
polycyclic  molecules.  It  is  shown  that  the  proposed  potential  model  is  able  to  accurately 
reproduce  the  crystallographic  structures  and  lattice  energies  (where  available)  of  these  crystals. 
For  the  majority  of  these  crystals,  the  best  agreement  with  experimental  structural  and  energetic 
data  is  obtained  in  those  cases  when  the  electrostatic  charges  have  been  determined  using  ab 
initio  methods  that  include  electron  correlations  effects  (namely  MP2  and  B3LYP).  The  use  of 
the  electrostatic  charges  calculated  at  the  Hartree-Fock  level  results  in  large  deviations  of  the 
predicted  lattice  energies  from  the  experimental  values.  The  deviations  of  the  lattice  energies  can 
be  significantly  decreased  by  scaling  the  electrostatic  charges  with  a  constant  factor. 
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1.  Introduction 


Atomistic  simulation  is  increasingly  gaining  acceptance  as  a  practical  research  tool  in  the 
investigation  of  the  behavior  of  condensed-phase  materials.  In  addition  to  providing  information 
that  is  difficult  or  impossible  to  measure,  prediction  leads  to  a  reduction  in  unnecessary 
measurement  or  synthesis  of  candidates  in  the  course  of  design  of  new  materials.  However,  the 
power  of  atomistic  simulation  can  only  be  realized  if  the  description  of  the  molecular  system  is 
accurate.  The  development  of  accurate  intermolecular  potentials  is  not  a  simple,  straightforward 
procedure.  Substantial  work  has  been  directed  toward  determining  both  simple  functions  that 
make  large-scale  simulation  realizable,  and  correct  parameterization  such  that  the  physical 
properties  of  the  materials  are  properly  described.  In  this  work,  an  intermolecular  potential  that 
accurately  describes  nitramine  crystals  is  presented.  Also,  the  potential  parameters  effect  on 
predictive  ability  is  investigated. 

In  initial  studies  of  nonreactive  processes  in  the  nitramine  explosive  hexahydro- 1,3,5, -trintro- 
1,3,5-s-triazine  (RDX),  an  intermolecular  potential  energy  function  was  developed  that  would 
accurately  reproduce  the  structure  of  the  a-form  of  the  RDX  crystal  [1].  This  potential  is 
composed  of  pairwise  atom — atom  (6-exp)  Buckingham  terms  with  explicit  inclusion  of  the 
electrostatic  interactions  between  the  charges  associated  with  the  atoms  of  different  molecules. 
The  parametrization  of  the  potential  function  was  done  such  that  molecular  packing  calculations 
(MP)  reproduced  the  experimental  structure  of  the  crystal  and  its  lattice  energy. 
Isothermal-isobaric  molecular  dynamics  simulations  (NPT-MD)  using  this  potential  energy 
function  predicted  crystal  structures  in  excellent  agreement  with  the  experimental  data  [1]. 

It  has  been  shown  that  this  interaction  potential  energy  function  is  transferable  to  two  other 
nitramine  crystals:  the  polycyclic  nitramine  2,4,6,8,10,12-hexanitrohexaazaisowurtzitane  (HNIW) 
[2]  and  the  monocyclic  nitramine  l,3,5,7-tetranitro-l,3,5,7-tetraazacyclooctane  (HMX)  [3].  Both 
MP  and  NPT-MD  simulations  predict  geometrical  parameters  in  good  agreement  with  the 
experimental  values  for  the  different  polymorphs  of  the  HNIW  and  HMX  crystals  [2,  3]. 
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Furthermore,  the  calculations  indicate  a  stability  ranking  for  HNIW  in  agreement  with 
experimental  measurements  [4]. 

The  success  of  these  potential  energy  parameters  in  describing  the  RDX  crystal  and  different 
phases  of  the  HMX  and  HNIW  crystals  provides  impetus  for  further  investigations  to  determine 
the  limits  of  the  transferability  of  this  interaction  potential.  Toward  this  end,  MP  calculations  of 
30  nitramine  crystals  are  reported  here.  This  set  of  crystals  is  composed  of  monocyclic, 
polycyclic,  and  acyclic  nitramine  molecules.  There  was  a  particular  interest  to  see  if  the 
geometrical  and  energetic  parameters  for  these  crystals  could  be  reproduced  by  this  proposed 
model. 

One  of  the  main  factors  that  contributes  to  the  quantitative  description  of  the  molecular 
packing  in  a  crystal  is  related  to  the  representation  of  the  electrostatic  interactions.  It  was  shown 
more  than  a  decade  ago  [5,  6]  that  increased  accuracy  in  structural  predictions  of  the  molecular 
crystals  and  in  transferability  of  the  potential  parameters  can  be  achieved  by  explicit  use  of  the 
electrostatic  interactions  between  the  charges  associated  with  the  atoms.  For  example,  many  of 
the  available  force  fields  such  as  Amber  [7],  ECEPP  [8],  or  Dreiding  [9]  use  these  kinds  of 
potential  terms  in  models  of  organic,  biological,  and  main  group  inorganic  crystals.  Further 
improvement  of  the  description  of  the  electrostatic  forces  between  molecules,  particularly  in 
crystals  with  substantial  anisotropies,  can  be  achieved  by  using  sets  of  point  multipoles  (charge, 
dipole,  quadrupole,  etc.)  on  every  atomic  site.  This  distributed-multipole  representation  has  been 
shown  to  be  successful  in  the  modeling  of  the  crystal  structures  of  polar  and  hydrogen-bonded 
molecules  [10]. 

The  present  study  has  found  that,  as  in  the  cases  of  the  RDX  [1],  HNIW  [2],  and  HMX  [3] 
crystals,  the  set  of  30  crystals  considered  here  can  be  accurately  represented  using  the 
Buckingham  potential  plus  Coulombic  interactions.  The  assignment  of  the  electrostatic  charges 
poses  a  problem  in  that  the  atom-centered  monopole  charge  is  not  an  observable  quantity  and 
cannot  be  obtained  directly  from  either  experiment  or  first  principles  calculations.  Currently, 
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there  are  several  schemes  for  evaluation  of  charges  by  empirical  partition  or  by  using  a  quantum 
mechanically  derived  wave  function  [11-13].  The  Coulombic  terms  have  been  determined 
through  fitting  of  partial  charges  centered  on  each  atom  of  the  molecules  to  a  quantum 
mechanically  derived  electrostatic  potential  [13].  An  investigation  has  been  done  on  how  the 
geometrical  and  energetic  parameters  predicted  in  MP  calculations  depend  on  charges  determined 
from  ab  initio  methods,  which  do  or  do  not  include  electron  correlation  effects.  Specifically, 
different  sets  of  charges  derived  from  the  Hartree-Fock  (HF)  wave  function  [14]  and  from 
methods  that  employ  electron  correlations  such  as  second-order  Moller-Plesset  (MP2) 
[15-18]  and  B3LYP  [19,  20]  have  been  used. 

The  studies  described  here  represent  the  first  stage  in  the  development  of  a  general  model  for 
nitramine  crystals.  The  main  limitations  of  the  present  model  are  due  to  the  assumption  of  rigid 
molecules  but  further  refinement  of  this  model  can  be  made  to  include  the  effects  of 
intramolecular  motions,  particularly  of  low-frequency  torsional  motions  of  the  nitro  groups  and 
the  ring. 

The  organization  of  the  paper  is  as  follows.  In  section  2,  the  intermolecular  potential  used  to 
simulate  the  nitramine  crystals  is  presented.  In  sections  3  and  4,  the  molecular  packing  methods 
and  results,  respectively,  are  described.  The  main  conclusions  are  summarized  in  section  5. 

2.  Intermolecular  Potential 


The  central  problem  in  classical  simulations  of  molecular  crystals  is  the  construction  of 
realistic  potentials  that  accurately  predict  the  structural  and  thermochemical  parameters.  In  this 
paper,  the  same  general  model  is  employed  for  the  atom-atom  potentials  that  proved  to  be 
successful  in  modeling  of  the  RDX,  HNIW,  and  HMX  crystals  [1-3].  In  particular,  it  is  assumed 
that  (1)  the  intermolecular  interactions  depend  only  on  the  interatomic  distances;  (2)  the 
interaction  potential  can  be  separated  in  contributions  identified  as  van  der  Waals  and 
electrostatic,  and  (3)  the  same  type  of  van  der  Waals  potential  is  used  for  the  same  type  of  atoms, 
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independent  of  their  valence  state.  Moreover,  in  the  present  case,  the  transferability  of  the 
potential  parameters  determined  for  the  case  of  RDX  crystal  to  all  the  nitramines  considered  in 
this  study  is  assumed. 

In  the  present  treatment,  the  intermolecular  interactions  between  the  molecules  of  the  crystal 
as  the  sum  of  pairwise  Buckingham  (6-exp)  (repulsion  and  dispersion)  and  Coulombic  (C) 
potentials  is  approximated: 


V 6-exp 
vap 


(r)=Aopexp(-Bapr)-Cap/r6, 


(1) 


and 


4718  0r  ’ 


(2) 


where  r  is  the  interatomic  distance  between  atoms  a  and  p,  qa  and  qp  are  the  electrostatic  charges 
on  the  atoms,  and  so  is  the  dielectric  permittivity  constant  of  vacuum. 


The  parameters  for  the  6-exp  potential  in  equation  (1)  are  those  previously  determined  for  the 
case  of  the  RDX  crystal  [1].  The  same  combination  rules  are  used  for  the  calculation  of  the 
heteroatom  parameters  from  homoatom  parameters,  as  previously  reported  [1], 

The  assignment  of  the  electrostatic  charges  is  made  by  using  the  set  of  atom-centered 
monopole  charges  for  the  isolated  molecule  that  best  reproduces  the  quantum  mechanically 
derived  electrostatic  potential,  which  is  calculated  over  grid  points  surrounding  the  van  der 
Waals  surface  of  the  molecules.  This  method  of  fitting  the  electrostatic  potential  was  proposed 
by  Breneman  and  Wiberg  [13]  and  is  incorporated  in  the  Gaussian  94  package  of  programs  [21] 
under  the  keyword  CHELPG  (electrostatic-potential-derived  atomic  charges).  This  method  has 
the  advantage  of  a  higher  density  of  points  and  a  better  selection  procedure,  which  ensures  a 
significant  decrease  in  orientation  effects  compared  to  those  observed  with  similar  methods  [12]. 
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The  CHELPG  charges  were  found  to  be  invariant  to  either  the  rotation  of  the  molecular 
coordinates  or  internal  bond  rotations.  These  calculations  have  been  done  at  both  the  HF  [14] 
and  second-order  MP2  [15-18]  levels  to  investigate  the  effect  of  electron  correlation. 

For  the  purpose  of  comparison  and  as  an  alternative  to  the  computationally  demanding  MP2 
method,  density  functional  theory  (DFT)  has  also  been  used  in  the  Kohn-Sham  formulation  [22], 
The  DFT  methods  offer  a  less  expensive  but  still  accurate  computational  alternative  to  ab  initio 
methods  for  including  the  electron  correlation  in  post-HF  treatments.  In  particular,  the  exchange 
functional,  described  by  the  fitted  three-parameter  hybrid  of  Becke  [16]  and  the  correlation 
functional  of  Lee,  Yang,  and  Parr  (B3LYP)  [17]  was  employed.  All  the  aformentioned  theoretical 
calculations  were  done  using  a  reasonable  quality  basis  set  (i.e.,  6-3 1G**  [split- valence  plus  d- 
type  and  p-type  polarization  functions])  [24]. 

It  has  been  shown  previously  [25,  26]  that  the  neglect  of  electron  correlation  in  self-consistent 
wave  functions  overestimates  the  electrostatic  interactions;  however,  this  is  mainly  a  scaling 
effect.  Cox  and  Williams  [25]  have  suggested  that  a  scaling  factor  of  0.9  can  be  used  to  improve 
agreement  between  the  calculated  and  experimental  values  of  the  dipole  moments  for  a  set  of 
eight  small  molecules.  The  same  factor  has  been  justified  in  a  study  of  the  electrostatic 
interactions  of  a  dipeptide  [26],  as  well  as  in  a  more  recent  work  related  to  the  role  of 
electrostatic  interactions  in  determining  the  crystal  structures  of  polar  organic  molecules  [10]. 
Such  an  electrostatic  model  has  been  employed  to  further  evaluate  the  effects  of  this  scaling 
procedure.  Specifically,  four  electrostatic  models  were  tested  for  each  of  the  30  crystals.  Two  of 
them  use  electron  correlation  methods  (namely  MP2  and  B3LYP);  the  third  one  uses  unsealed  HF 
charges,  and  the  last  HF  charges  scaled  by  0.9  (denoted  as  0.9HF). 

3.  Computational  Approach 

A  general  procedure  for  testing  intermolecular  potential  energy  functions  for  organic  crystals 
is  based  on  the  use  of  molecular  packing  calculations  [5,  6].  The  basic  idea  is  to  minimize  the 
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lattice  energy  with  respect  to  the  structural  degrees  of  freedom  of  the  crystal.  For  crystals  with 
one  molecule  in  the  asymmetric  unit  occupying  an  arbitrary  position,  the  maximum  number  of 
degrees  of  freedom  is  12  and  corresponds  to  the  six  unit  cell  constants  (a,  b,  c,  a,  3,  y),  the  three 
rotations  (0i,  02,  03)  and  the  three  translations  (xi,  x2,  x3)  of  the  rigid  molecule.  A  reduced 
number  of  structural  degrees  of  freedom  might  be  involved,  depending  on  the  symmetry 
restrictions  of  different  space  groups.  For  crystals  with  more  than  one  molecule  in  the 
asymmetric  unit,  additional  degrees  of  freedom  are  introduced  to  describe  the  rotation  and 
translation  of  the  additional  molecules. 

Assuming  that  the  crystal  energy  is  known  as  a  function  of  the  structural  lattice  parameters, 
the  equilibrium  crystal  configuration  is  determined  by  the  conditions  of  zero  force  and  torque, 
together  with  the  requirement  that  there  is  a  minimum.  The  search  for  such  a  minimum  can  be 
done  using  a  combination  of  steepest-descent  and  Newton-Raphson  procedures  [27, 28]. 

In  the  present  study,  it  is  assumed  that  the  crystals  can  be  represented  as  an  ensemble  of  rigid 
molecules.  The  minimization  of  the  lattice  energy  with  symmetry  constraints  has  been 
performed  using  the  molecular  packing  program  PCK91  [29]  by  taking  the  experimentally 
observed  geometries  as  starting  configurations.  This  program  employs  an  accelerated 
convergence  method  [1,  28]  for  accurate  evaluation  of  the  crystal  Coulombic  and  dispersion 
lattice  sums,  with  the  first  and  second  derivatives  of  the  crystal  energy  evaluated  analytically.  In 
all  calculations,  a  cutoff  distance  of  19  A  has  been  used  with  the  parameters  T|  that  determines 
the  relative  contributions  of  the  real-  and  reciprocal-space  terms  as  defined  in  Sorescu  et  al.  [1] 
having  the  values  t|i  =  0.1861  A  1  and  T|6=  0.2304  A  L  The  space-group  symmetry  is  maintained 
throughout  the  energy  minimization.  This  reduces  the  number  of  independent  variables  in  the 
minimization  procedure,  resulting  in  a  significant  decrease  in  the  computational  time  compared 
to  unconstrained  energy  minimization.  For  example,  for  the  3-phase  of  HMX  crystal  with  space 
group  P2]/n  (Z  =  2),  the  crystallographic  parameters  that  were  varied  in  the  minimization  using 
the  PCK91  program  are  the  three  dimensions  of  the  unit  cell  and  the  angle  3,  while  the  angles 
ct  and  y  were  frozen  at  90°.  Since  the  molecule  in  the  asymmetric  unit  occupies  an  inversion 
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center,  only  the  three  rotations  of  this  molecule  were  allowed  to  vary,  while  the  three  translations 
were  not  modified  to  maintain  the  symmetry  imposed  by  the  inversion  center.  It  has  been 
previously  shown  [1,  2,  30]  that,  despite  the  symmetry  restrictions  imposed  in  the  PCK91 
program,  the  final  lattice  energies  and  crystallographic  parameters  are  in  good  agreement  with 
those  obtained  when  symmetry  constraints  were  removed  (i.e.,  the  predicted  crystals  maintained 
the  observed  space-group  symmetry). 

The  quality  of  the  predicted  geometrical  crystallographic  parameters  relative  to  the 
experimental  values  has  been  done  using  a  structural  shift  factor  of  the  form  [28, 31] 


F  =  (A0/2)2  +  (10Ax)2  +  (100Aa/a)2  +  (100Ab/b)2  +  (100Ac/c)2  +  (Act)2  +  (A0)2  +  (Ay)2,  (3) 

where  A0  is  the  total  root-mean-square  (rms)  rigid-body  rotational  displacement  (in  degrees)  after 
minimization;  Ax  is  the  rms  total  rigid-body  translational  displacement  (in  angstroms);  and  a,  b,  c, 
and  a,  p,  y  are,  respectively,  the  lengths  of  the  edges  and  the  angles  of  the  unit  cell. 

An  important  test  of  the  validity  of  the  model  is  the  accuracy  of  the  predicted  lattice  energies 
of  the  crystals.  The  lattice  energies  determine  the  relative  stabilities  of  the  different 
crystallographic  phases.  The  calculated  static  lattice  energy  can  be  compared  to  the  experimental 
sublimation  enthalpy  by  the  using  the  relationship  [32]  -  AHsubi  =  E  +  Ko  +  2RT,  where  E  is  the 
lattice  energy  and  Ko  is  the  zero-point  energy.  Often,  a  rough  estimation  of  the  lattice  energy  is 
obtained  by  neglecting  the  Ko  term.  Kitaigorodski  [5]  has  pointed^  out  that,  considering  the 
inaccuracy  involved  in  the  experimental  determination  of  AHa*  and  with  the  neglect  of  zero-point 
energy,  discrepancies  up  to  3-4  kcal/mol  between  the  calculated  and  the  observed  enthalpies  of 
sublimation  are  expected  [5].  In  the  case  of  RDX  crystals,  it  was  found,  using  this  approximation, 
that  the  predicted  lattice  energy  (E  =  -130.09  kJ/mol)  is  in  very  close  agreement  with  the 
experimental  sublimation  enthalphy  (-AHsub  =  - 130.1  kJ/mol)  [1]. 
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4.  Results  and  Discussions 


The  30  nitramine  molecules  considered  in  this  study  are  shown  in  Figure  1.  They  were  chosen 
as  representative  examples  of  important  acyclic  and  cyclic  nitramines.  Different  types  of  mono- 
and  polycyclic  nitramines  have  been  included,  particularly  crystals  that  are  important  energetic 
materials.  The  structures  of  most  of  these  crystals  have  been  obtained  by  x-ray  diffraction 
techniques.  Despite  the  generally  poorer  resolution  of  hydrogen  atom  positions  obtained  by  these 
techniques,  no  additional  adjustment  of  these  positions  has  been  done  to  give,  for  example,  the 
standard  bond  lengths  [33].  The  crystal  structures  in  Figure  1  are  denoted  first  by  the  common 
names  of  the  molecules.  Where  available,  the  crystal  abbreviation  from  the  original  reference  is 
also  included  with  the  corresponding  crystal  “refcode”  used  in  the  Cambridge  Structural  Database 

[34]  and  indicated  by  the  term  in  the  second  set  of  parentheses.  The  structures  used  for  HNIW 

[35]  and  fi-HMX  [36]  crystals  are  not  in  the  Cambridge  database,  so  they  do  not  have  a  refcode. 
In  addition,  different  crystallographic  phases  of  HMX  and  HNIW  crystals  have  been  studied; 
these  are  detailed  in  Table  1.  The  specific  references  for  each  of  the  30  crystals  in  Table  1  are 
provided  in  brackets  in  the  “Crystal”  column  and  found  in  the  reference  section  [35-57]. 

The  results  of  MP  calculations  using  the  PCK91  program  are  presented  in  Table  1.  The 
predicted  structural  lattice  parameters  for  the  great  majority  of  these  crystals  deviate  by  less  than 
2%  from  the  experimental  structures.  Also,  for  the  majority  of  the  crystals,  there  are  small 
rotations  and  practically  no  translations  for  the  molecules  in  the  asymmetric  unit  cell.  The 
accuracy  of  the  predictions  can  be  seen  in  Figure  2,  where  the  overall  structural  drift  factors 
described  in  equation  (3)  are  given.  Only  10%  of  the  total  number  of  crystals  considered  here 
have  a  structural  shift  factor  larger  than  2.0,  and  practically  half  of  the  crystals  have  shift  factors 
that  are  less  than  1.0. 

It  is  important  to  point  out  that,  ideally,  the  predicted  lattice  structural  parameters  should  be 
compared  with  the  values  determined  at  zero  temperature.  However,  this  is  not  possible  due  to 
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Figure  1.  Illustration  of  Molecules  Whose  Crystal  Structures  Were  Studied.  The  Common 
Abbreviation  of  the  Crystal  Name  Is  Given  in  the  First  Set  of  Parentheses,  and 
the  Refcode  Entry  of  the  Cambridge  Structural  Database  [34]  Is  Given  in  the 
Second  Set  of  Parentheses. 
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Table  1.  Comparison  of  the  Crystallographic  Parameters  as  Determined  in  Molecular  Packing  Calculations  to  the 
Corresponding  Experimental  Values  for  Different  Molecular  Structures 


,  values  in  parentheses  represent  the  procentual  difference  relative  to  experimental  values. 

For  monoclinic  systems,  only  angles  are  given.  For  triclinic  systems,  the  cell  angles  and  the  cell  angle  errors  are  given  in  the  order.  For  orthorhombic  crystal 
systems,  a  value  of  90.0  is  indicated  only  for  the  experimental  structure. 

Jx  =  the  net  rotation  of  the  molecule  and  the  net  translation  of  the  center  of  mass  of  the  molecule. 

F  =  structural  shift  factor;  defined  in  equation  (3). 
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Table  1.  Comparison  of  the  Crystallographic  Parameters  as  Determined  in  Molecular  Packing  Calculations  to  the 
Corresponding  Experimental  Values  for  Different  Molecular  Structures  (continued) 


i  values  in  parentheses  represent  the  procentual  difference  relative  to  experimental  values. 

For  monoclinic  systems,  only  angles  are  given.  For  triclinic  systems,  the  cell  angles  and  the  cell  angle  errors  are  given  in  the  order.  For  orthorhombic  crystal 
systems,  a  value  of  90.0  is  indicated  only  for  the  experimental  structure. 

5x  =  the  net  rotation  of  the  molecule  and  the  net  translation  of  the  center  of  mass  of  the  molecule. 

F  =  structural  shift  factor;  defined  in  equation  (3). 
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Table  1.  Comparison  of  the  Crystallographic  Parameters  as  Determined  in  Molecular  Packing  Calculations  to  the 
Corresponding  Experimental  Values  for  Different  Molecular  Structures  (continued) 


a  values  in  parentheses  represent  the  procentual  difference  relative  to  experimental  values. 

For  monoclinic  systems,  only  angles  are  given.  For  triclinic  systems,  the  cell  angles  and  the  cell  angle  errors  are  given  in  the  order.  For  orthorhombic  crystal 
systems,  a  value  of  90.0  is  indicated  only  for  the  experimental  structure. 
bx  =  the  net  rotation  of  the  molecule  and  the  net  translation  of  the  center  of  mass  of  the  molecule. 

CF  =  structural  shift  factor;  defined  in  equation  (3). 
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For  monoclinic  systems,  only  angles  are  given.  For  triclinic  systems,  the  cell  angles  and  the  cell  angle  errors  are  given  in  the  order.  For  orthorhombic  crystal 
systems,  a  value  of  90.0  is  indicated  only  for  the  experimental  structure. 

bx  =  the  net  rotation  of  the  molecule  and  the  net  translation  of  the  center  of  mass  of  the  molecule. 

CF  =  structural  shift  factor;  defined  in  equation  (3). 
dFor  P31  and  P61,  the  angle  =  120  is  indicated. 
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Crystal  Index 

Figure  2.  The  Calculated  Structural  Shift  Factor  F  (Equation  [31)  for  the  Crystal 
Structures.  The  Crystal  Index  Corresponds  to  the  Number  of  the  Crystal  in 
Table  1.  The  Horizontal  Lines  at  1  and  2  Are  Marked  for  a  Clearer  View  of  the 
Distribution  of  Points. 

lack  of  data  at  low  temperatures.  Consequently,  the  aforementioned  comparison  considers  the 
deviations  of  the  predicted  geometrical  parameters  from  the  experimental  values  obtained  at 
room  temperature.  It  can  be  observed  from  the  data  given  in  Table  1  that  the  predicted  lattice 
dimensions  either  underestimate  or  overestimate  the  experimental  values.  Consequently,  there  is 
not  a  general  trend  of  the  relationship  between  the  predicted  and  the  experimental  geometric 
lattice  periods,  despite  the  small  deviations  between  the  two  sets  of  values. 

The  influence  of  the  level  of  ab  initio  calculations  on  the  final  crystallographic  parameters  is 
also  illustrated  by  the  results  in  Table  1.  The  difference  in  the  predicted  geometrical  parameters 
is  less  than  1 .0%  when  correlated  and  uncorrelated  ab  initio  methods  are  used.  In  addition,  there 
is  not  a  clear  trend  of  the  degree  of  accuracy  with  the  ab  initio  level  of  calculations.  In  1 9  of  the 
systems,  the  accuracy  increases  when  the  HF  charges  are  replaced  with  those  obtained  at  the 
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B3LYP  and  MP2  levels,  but,  in  the  other  1 1  cases,  the  accuracy  decreases.  The  maximum 
difference  of  the  structural  shift  factors  with  the  different  levels  of  calculation  is  less  than  0.27%. 

It  can  also  be  seen  from  the  results  in  Table  1  that,  when  the  electrostatic  charges  calculated 
at  HF  level  are  scaled  by  the  0.9  factor,  the  corresponding  predicted  geometrical  parameters  are 
very  close  to  those  obtained  at  MP2  level.  Moreover,  the  corresponding  structural  shift  factors 
have  values  intermediate  between  the  MP2  and  HF  values. 

The  lattice  energies  predicted  by  different  models  are  given  in  Table  2.  As  can  be  seen  by 
comparing  the  results  for  MP2,  B3LYP,  and  HF  methods,  the  use  of  the  correlated  methods 
result  in  a  decrease  in  the  absolute  lattice  energy.  This  effect  can  be  understood  as  a  consequence 
of  the  decrease  in  the  absolute  value  of  the  electrostatic  energy,  which  is  attractive.  The 
variations  in  the  absolutes  values  of  the  HF  lattice  energies  are  between  8.5  and  17.5%,  relative 
to  the  MP2  energies,  with  the  average  deviation  12.8%  (see  Figure  3).  The  use  of  the  0.9  scaling 
factor  reduces  these  deviations  to  the  range  0-7.8%,  with  the  average  deviation  4.1%.  Finally, 
the  B3LYP  lattice  energies  are,  as  expected,  much  closer  to  the  MP2  energies,  with  the  range  of 
variations  1.5-3 .9%  and  average  deviation  2.6%.  These  results  indicate  that  the  lattice  energies 
differ  significantly  for  sets  of  electrostatic  charges  calculated  with  ab  initio  methods  that  do  not 
include  electron  correlation.  These  differences  can  be  decreased  by  a  factor  of  ~3  by  scaling  the 
HF  charges.  Another  important  result  is  that  DFT  can  provide  charges  that  give  an  accuracy 
(within  2.6%)  for  the  lattice  energy  that  is  comparable  to  those  determined  at  the  MP2  level. 
These  results  are  important  since  the  computational  times  for  B3LYP  calculations  are 
significantly  lower  than  those  for  MP2. 

The  calculated  lattice  energies  can  also  be  compared  to  experimental  sublimation  enthalpies 
in  Table  2.  For  RDX  (CTMTNA),  HMX  (0-HMX,  OCHTET03)  and  dimethylnitramine 
(METNAM08)  crystals  the  agreement  of  the  MP2  energies  to  the  experimental  values  is  very 
good,  while  for  tetryl  (NTNANL)  the  difference  of  15.7  kJ/mol  is  within  the  range  12-17  kJ/mol 
considered  acceptable  by  Kitaigorodsky  [5].  Despite  the  limited  number  of  experimental  values 
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Table  2.  Comparison  of  the  Experimental  and  Calculated  Lattice  Energies  for  Different 
Sets  of  Electrostatic  Charges 


IncT 

Crystal 

AHsub 

(kJ/mol) 

Lattice  Energies 
(kJ/mol) 

MP2/6-31G** 

B3LYP/6-31G** 

0.9  HF/6-31G** 

HF/6-31G** 

T 

CTMTNA 

130.1  [58] 

-130.09 

-133.68 

-136.05 

-148.18 

2 

CIWMEA10 

-113.01 

-116.40 

-118.81 

-129.07 

3 

KOFKAR 

-116.55 

-119.40 

-119.61 

-129.26 

4 

NOHTAZ 

-141.45 

-145.35 

-145.99 

-157.25 

5 

KOFKEV 

-120.00 

-122.66 

-122.20 

-131.95 

6 

p-HMX 

175.2  [58] 

-180.23 

-185.02 

- 187.01 

-204.45 

7 

OCHTET 

-179.15 

-184.66 

-188.74 

-208.75 

8 

OCHTET03 

161.9  [59] 

-168.24 

-173.44 

-178.11 

-197.67 

9 

CATJIQ 

-182.50 

-187.59 

-196.88 

-213.26 

10 

JEXLUT 

-178.54 

-184.32 

-186.70 

-200.43 

11 

JEXMH 

-169.22 

-172.11 

-177.36 

-190.98 

12 

JEXMEE 

-172.98 

-178.48 

-179.86 

-196.13 

13 

KEMTIF 

-184.42 

- 190.46 

-192.56 

-208.80 

14 

JEXMAA 

-176.33 

-182.29 

181.97 

-196.64 

15 

SECVOL 

-187.66 

- 192.30 

-193.63 

-209.42 

16 

HNIW  (e-phase) 

-186.77 

- 192.82 

-196.52 

-210.88 

17 

HNIW  (p-phase) 

-181.29 

- 186.28 

-190.54 

-201.81 

18 

HNIW  (y-phase) 

-175.31 

- 180.90 

-186.15 

-198.61 

19 

DNPMTA 

-140.60 

-142.78 

-142.60 

- 154.30 

20 

MTNANL 

133.8  ±  1.6  [60] 

-149.53 

-153.81 

-158.30 

-170.57 

21 

KOFKIZ 

-262.18 

-267.98 

-265.78 

-288.98 

22 

JEDSUG 

-124.22 

-128.89 

-131.24 

-145.14 

23 

JEHLAJ 

-180.70 

-187.76 

-192.29 

-209.63 

24 

METNAM08 

69.87  [61] 

-70.26 

-71.20 

-70.21 

-76.25 

25 

GEJXAU 

-171.87 

-174.31 

-175.69 

-191.41 

26 

NXENAM01 

-131.84 

-134.80 

-137.57 

-148.40 

27 

NABMUY01 

-169.28 

-173.70 

-180.86 

- 196.26 

28 

DILFUZ 

-216.09 

-220.50 

-209.57 

-241.50 

29 

NOETNA02 

-166.55 

-163.68 

-172.29 

-181.04 

30 

ACENIH 

-186.20 

-188.32 

-191.61 

-205.65 

available  for  comparisons,  it  can  be  seen  that  a  significant  improvement  in  the  accuracy  of  the 
predicted  lattice  energies  can  be  obtained  by  using  the  electrostatic  charges  determined  by 
electron  correlated  methods.  The  scaling  of  the  HF  charges  also  leads  to  improvements  of  the 
predicted  energies,  but  the  differences  from  the  experimental  values  are  larger  than  those  obtained 
when  the  charges  are  calculated  with  electron  correlation  methods. 
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Figure  3.  The  Percent  Differences  Between  the  Lattice  Energies  and  Those  Based  on  the 
MP2  Values.  The  Crystal  Index  Corresponds  to  the  Number  of  the  Crystals  in 
Table  1.  The  Three  Horizontal  Lines  Indicate  the  Average  Deviations  for  the 
Energies  Calculated  Using  the  B3LYP  (<pl>  =  2.6%),  0.9*HF  (<p2>  =  4.1%) 
and  HF  (<p3>  =  12.8%)  Sets  of  Charges. 
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The  relative  stability  of  different  polymorphic  phases  of  the  HMX  and  HNIW  crystals  has  also 
been  investigated.  The  calculated  MP2  lattice  energies  for  the  0,  a,  and  8  phases  of  HMX  are 
-180.23,  -179.15  and  -168.24  kJ/mol,  respectively.  These  values  support  the  polymorph 
stability  ranking  (3  >  a  >  8  found  experimentally  by  McCrone  [62].  Also,  the  calculated  lattice 
energies  per  molecule  for  the  e-,  0-,  and  y-HNIW  phases  of  -186.77,  -181.29,  and 
-175.31  kJ/mol,  respectively,  are  consistent  with  the  stability  ranking  e  >  0  >  y  reported  by 
Russell  et  al.  [4], 

5.  Conclusions 

An  investigation  has  been  done  on  the  transferability  of  a  6-exp  Buckingham  potential 
previously  developed  for  the  a-RDX  crystal  [1]  to  30  crystals,  consisting  of  acyclic,  monocyclic, 
or  polycyclic  nitramines.  The  intermolecular  potential  includes  Coulombic  interactions  between 
electrostatic  charges.  These  charges  were  determined  from  fits  to  ab  initio  electrostatic  potentials 
calculated  for  the  individual  molecules  in  the  experimental  configurations. 

The  tests  of  this  potential  for  the  set  of  30  crystals  have  been  performed  using  molecular 
packing  calculations.  Accurate  values  of  the  crystal  lattice  energy  have  been  obtained  by 
employing  the  accelerated  convergence  technique  for  the  dispersion  and  Coulombic  lattice  sums. 
Four  different  electrostatic  models,  have  been  considered  with  charges  determined  at  HF,  B3LYP, 
and  MP2  levels,  and  with  charges  obtained  at  HF  level  uniformly  scaled  by  a  factor  of  0.9.  The 
predicted  geometries  indicate  a  good  agreement  with  the  experimental  values  for  the  great 
majority  of  the  crystals  in  the  study.  For  90%  of  the  crystals,  the  structural  shift  factor  was  less 
than  2.0,  while  for  50%  of  them,  it  is  less  than  1.0. 

There  is  only  a  small  influence,  generally  below  1%,  on  the  crystallographic  parameters  by  the 
set  of  electrostatic  charges  used.  However,  the  lattice  energies  are  strongly  dependent  on  the 
electrostatic  model.  In  particular,  the  best  overall  agreement  with  the  experimental  lattice 
energies  was  obtained  by  using  MP2  calculated  charges.  The  lattice  energies  calculated  using  the 
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B3LYP  charges  overestimate  the  MP2  energies  by  about  2.6%,  while  the  overestimation  in  the 
case  of  HF  charges  is  about  12.8%.  The  procedure  of  uniformly  scaling  the  HF  charges,  [10, 25] 
decreases  the  differences  to  about  4.1%.  It  was  also  shown  that  this  intermolecular  potential 
correctly  describes  the  order  of  stability  of  different  phases.  The  predicted  stability  P  >  a  >  5  is 
in  accord  with  the  experimental  findings  [62].  Also,  the  calculated  stability  ranking  £  >  P  >  y  for 
HNIW  agrees  with  the  previously  reported  results  [4]. 

The  success  of  the  present  potential  energy  parameters  in  describing  different  types  of 
nitramines  and  different  phases  at  moderate  temperatures  and  low  pressure  provides  incentive  to 
further  investigate  the  transferability  of  this  model  to  other  classes  of  crystals.  Incorporation  of 
intramolecular  motion  by  relaxing  the  rigid  molecular  model  will  also  be  investigated  in  future 
studies. 
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